% Done late at night in 2021,as I reanalyze 4 year old data.
% Load data

%% Determine
clear all
fclose all
path=['C:\Users\Tom�s\Desktop\Temp\'];

fileNumberNoDrugs=4; %First file sigh. what a stupid way to do this, but I was young.
trials=3;



ISteps=[-30:5:30];

lookup=[7,8,9,10,11,12,13,6,5,4,3,2,1]; %%2021--> This is confusing AF, but I got it. The number indicates the index of the corresponding value of current injection in the vector -30:5:30. e.g. step_7 -->lookup(7+1)=6--> ISteps(6)=-5 pA  Therefore step_7 corresnpoinds to -5 pA. fk lol

tracesNoDrugs=NaN(550000,numel(ISteps),trials);



    for i=1:numel(ISteps)
        for j=1:trials;
       Ibw1= IBWread([path 'ch0_' num2str((fileNumberNoDrugs-1)+i+(13*(j-1))) '.ibw']);
      
       idxND1=strfind(Ibw1.WaveNotes,'step_');
        idxND2=strfind(Ibw1.WaveNotes,';HC');
      
        tracesNoDrugs(:,lookup(str2num(Ibw1.WaveNotes(idxND1+5:idxND2-1))+1),j)=medfilt1(Ibw1.y,10);
        
    end 
    end


derivative=diff(tracesNoDrugs,2);
derivative=mean(derivative,3); %average all trials


meandiff=median(abs(derivative),2);
diffMAD=mad(meandiff(1:1e3));
diffBaseline=mean(meandiff(1:1e3));
threshCrossing=meandiff>(diffBaseline+(20*diffMAD));
StartStepIdx = find(threshCrossing, 1, 'first'); %find index of first crossing of 10 MADs of mean derivative of steps. This corresponds with the start of the step.
StartStepIdx=1e5; %override for when it's whacky and I can just set it manually
EndIdx=StartStepIdx+ (0.25/Ibw1.dx)% Index of 250ms after step start. This is the region for counting spikes for IO curves regardless of step duration.



Dx=Ibw1.dx;


for i=1:numel(ISteps)
    for j=1:trials
    
    %find spikes and count them.
    diffTemp=medfilt1(diff(tracesNoDrugs(:,i,j)),10);
    diffMADtemp=mad(diffTemp(1:1e3));
    [pks,locs]=findpeaks(tracesNoDrugs(:,i,j),'minPeakprominence',15,'minpeakDistance',200);
    
    numSpikes(i,j)=numel(locs(locs<EndIdx & locs>(StartStepIdx+100))); % only count first 250 ms
    
    FRSteps(i,j)=1/(Dx*(mean(diff(locs(locs<EndIdx)))));
    
    figure;plot(tracesNoDrugs(:,i,j))
xlim([StartStepIdx-10e3 EndIdx+ 10e3])
hold on
plot(locs,pks,'ro')
vline(StartStepIdx,'--k')
vline(EndIdx,'--k')

    
    
    
    
    
    
    end
end


FRSteps(isnan(FRSteps))=0; %get rid of the nans for plotting

figure;plot(ISteps ,numSpikes,'linewidth',2)
hold on 
plot(ISteps,numSpikes,'k.','markersize',10)
ylabel('spike number')
xlabel('current injection (pA)')

% figure;plot(Step,FRSteps)
% hold on 
% plot(Step,FRSteps,'ko')
% ylabel('firing rate (Hz)')
% xlabel('current injection (pA)')
% 
% figure;plot(tracesNoDrugs(:,i,j))
% xlim([StartStepIdx-10e3 EndIdx+ 10e3])
% hold on
% plot(locs,pks,'ro')
% vline(StartStepIdx,'--k')
% vline(EndIdx,'--k')

meanSpikes=(mean(numSpikes,2))'


fclose all